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1.  SUMMARY 


The  micromechanical  damage  mechanics  formulated  by  Ashby  and  Sainmis  (1990)  was  used  to 
explore  the  effect  of  ice  in  the  cracks  of  crystalline  rock  on  the  seismic  coupling  of  nuclear 
explosions  detonated  in  the  permafrost  layer  at  arctic  test  sites.  Experimental  data  show  that 
frozen  rock  is  elastically  stiffer  and  has  a  higher  fracture  strength  that  either  saturated  or  dry  rock 
above  the  freezing  point  of  water,  and  that  these  increases  in  stiffness  and  strength  increase  as 
temperature  is  decreased.  The  damage  mechanics  gave  a  good  description  of  these  data  if  the 
principal  effect  of  ice  in  the  cracks  is  to  increase  the  coefficient  of  sliding  friction.  The  known 
temperature  dependence  of  creep  in  ice  gave  a  quantitative  description  of  the  increase  in  strength 
with  decreasing  temperature.  An  important  implication  of  this  model  is  that  the  strength  should 
also  be  strain-rate  dependent,  making  it  an  important  factor  even  at  temperatures  near  the  melting 
point  of  ice  due  to  the  high  loading  rates  in  the  compressive  front  of  an  explosion. 

The  damage  mechanics  model  was  also  used  to  help  interpret  the  results  of  a  field  test  conducted 
in  collaboration  with  Weston  Geophysical  and  New  England  Research.  Identical  chemical 
explosions  were  detonated  in  frozen  and  unfrozen  schist  near  Fairbanks  Alaska.  Explosions  at 
the  frozen  site  produced  more  high-frequency  energy  resulting  in  a  higher  corner-frequency  in 
the  seismic  spectrum  and,  by  implication,  a  smaller  elastic  radius.  This  result  is  consistent  with 
damage  mechanics  simulations  in  which  ice  in  the  cracks  suppresses  fracture  damage  and 
predicts  a  smaller  elastic  radius.  In  these  simulations,  ultrasonic  velocity  measurements  by  New 
England  Research  on  frozen  and  unfrozen  lab  samples  from  the  two  sites  were  used,  together 
with  seismic  velocities,  to  calculate  the  insitu  initial  damage  at  the  frozen  and  unfrozen  test  sites. 


2.  INTRODUCTION 

The  overall  objective  of  this  research  program  has  been  to  use  micromechanical  damage 
mechanics  developed  by  Ashby  and  Sainmis  (1990)  to  describe  the  effects  of  the  extensive 
fracture  and  granulation  of  rock  known  to  occur  in  the  non-linear  source  region  of  underground 
nuclear  explosions.  The  initial  objective  was  to  explain  the  anomalously  wide  seismic  pulses 
produced  by  explosions  in  crystalline  rock  that  were  not  generated  by  the  numerical  simulations. 
The  introduction  of  fracture  damage  into  the  numerical  models  has  removed  this  discrepancy 
between  theory  and  the  observation,  and  the  micromechanical  model-based  damage  rheology 
developed  by  Ashby  and  Sainmis  (1990)  has  been  built  into  the  finite  difference  computer 
models  currently  used  to  model  the  seismic  coupling  of  underground  nuclear  explosions  (Rimer 
et  ah,  1999;  Stevens  et  ah,  2002). 

One  interesting  consequence  of  having  a  micromechanical  model  for  the  damage  is  that  the 
seismic  radiation  from  the  individual  fractures  can  be  calculated.  Although  the  radiation  from 
each  fracture  is  small  and  very  high-frequency,  Johnson  and  Sainmis  (2001)  found  that  the  net 
effect  of  integrating  all  the  fractures  is  the  radiation  of  significant  secondary  P  wave  energy  in 
the  seismic  band.  If  there  is  a  preferred  regional  pre-stress  direction,  or  if  the  initial  damage  has  a 
preferential  orientation,  or  both,  then  the  damage  process  can  also  radiate  significant  secondary  S 
wave  energy  in  the  seismic  band.  This  result  is  relevant  to  the  monitoring  program  because 
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recent  advances  in  seismic  discrimination  and  yield  estimate  of  underground  nuclear  explosions 
have  been  based  largely  on  high-frequency  local  phases  such  as  Lg  and  higher-mode  surface 
waves.  This  shift  in  focus  to  higher  frequencies  has  stimulated  new  interest  in  understanding  the 
non-linear  seismic  coupling  near  the  source.  One  result  of  our  current  phase  of  research  reported 
below  is  the  direct  observation  of  S  wave  generation  by  damage  in  high-speed  digital 
photographs  of  laboratory  explosions  in  photoelastic  plates. 

The  major  advantage  of  using  a  micromechanical  damage  model  is  that  the  effects  on  seismic 
coupling  of  depth  of  burial,  ground  water,  and  fracture  patterns  in  the  source  rock  can  be 
accounted  for  in  a  physically  meaningful  way  by  modeling  their  effect  on  the  nucleation,  growth 
and  interaction  of  fractures  in  the  source  rock.  For  example,  the  question  of  the  effect  of 
permafrost  on  seismic  coupling  has  recently  arisen.  The  bulk  of  this  report  details  our  use  of  the 
damage  mechanics  model  to  explore  the  effect  of  ice  in  cracks  on  seismic  coupling  in  frozen 
crystalline  rock.  We  first  use  the  model  to  explain  lab  data  on  the  effects  of  ice  on  strength  and 
elastic  stiffness  of  crystalline  rock.  We  then  model  the  results  of  a  field  experiment  carried  out  in 
collaboration  with  Weston  Geophysical  and  New  England  Research  in  which  identical  chemical 
explosions  were  detonated  in  frozen  and  un-frozen  gneiss  near  Fairbanks  Alaska. 

3.  TECHNICAL  APPROACH 

3.1  Non-linear  Regimes  in  the  Source  Region  of  an  Underground  Explosion 

Several  different  processes  occur  in  the  rock  near  a  tamped  underground  explosion  as  illustrated 
in  Figure  1  (see,  e.g.,  Rodean,  1971).  The  explosion  is  initially  contained  within  a  cavity  of 
radius  rc,  which  has  been  excavated  from  the  surrounding  rock.  At  the  time  of  detonation  a  hot 
pressurized  gas  is  created  within  the  cavity,  which  causes  it  to  expand.  Some  of  the  surrounding 
rock  may  be  vaporized  and  added  to  the  cavity  gas  at  this  time.  The  sudden  expansion  of  the 
cavity  generates  a  shock  wave  that  propagates  outward  causing  major  damage  to  the  surrounding 
rock  and  producing  the  series  of  regimes  summarized  in  Figure  1.  Rock  first  flows  plastically; 
then  is  stressed  beyond  its  brittle  failure  limit  and  becomes  granulated;  then  is  stressed  to  the 
point  where  radial  cracks  grow  but  failure  is  not  reached;  then  defonns  anelastically  where  pre¬ 
existing  cracks  slide  but  do  not  grow;  and  finally  deforms  elastically.  These  regimes  reflect  the 
decrease  in  energy  density  in  the  shock  with  distance  from  the  explosion  caused  partly  by 
spherical  spreading  and  partly  by  the  fact  that  energy  is  being  used  to  fracture  and  deform  the 
rock.  The  shock  wave  gradually  decays  into  an  inelastic  wave  involving  non-linear  motions, 
which  further  decays  with  distance  until  a  radius  is  reached  where  the  motions  are  small  enough 
to  be  described  by  the  ordinary  elastodynamic  equations  of  linear  elasticity.  Beyond  this  elastic 
radius,  re,  the  disturbance  caused  by  the  explosion  can  be  modeled  as  linear  elastic  waves  that 
propagate  throughout  the  rest  of  the  earth. 
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Figure  1.  Schematic  diagram  of  the  region  surrounding  a  contained  explosion  in  rock:  rc  is  the 
cavity  radius,  ry  is  the  radius  to  which  failure  occurs,  rj  is  the  radius  to  which  new 
damage  is  created,  and  re  is  the  elastic  radius  beyond  which  waves  can  be  approximated 
as  elastic. 

These  processes  that  occur  around  an  explosion  can  have  strong  effects  on  the  elastic  waves  that 
are  radiated  beyond  the  elastic  radius  re.  In  this  report  we  concentrate  on  the  growth  of  pre¬ 
existing  cracks  that  can  affect  the  radiated  elastic  waves  in  at  least  three  different  ways.  First, 
intense  cracking  will  significantly  lower  the  bulk  modulus  and  shear  modulus  near  the  source 
(O’Connell  and  Budianski,  1974;  Rimer  et  al.  1998).  Second,  it  has  been  speculated  that  when 
cracking  extends  into  the  failure  regime,  acoustic  fluidization  can  lower  the  basic  strength  of  the 
rock  (Sammis,  1998).  Third,  motions  on  the  cracks  serve  as  secondary  sources  of  elastic  waves 
that  can  contribute  to  the  net  seismic  radiation  field. 

The  spatial  extent  of  each  regime  and  the  associated  mechanics  are  calculated  using  the 
formulation  of  Ashby  and  Sammis  (1990)  to  model  the  nucleation  and  growth  of  cracks  in  the 
source  region,  which  is  here  referred  to  as  damage.  Calculation  of  this  damage  requires 
knowledge  of  the  stress  field  in  the  region  surrounding  the  source,  which  is  approximated  in  this 
study  by  the  equivalent  elastic  method  (Johnson,  1993;  Johnson  and  Sammis,  2001). 

3.2  Micromechanical  Damage  Mechanics 

The  concept  of  damage  used  here  is  that  developed  by  Ashby  and  Sammis  (1990).  In  that  paper, 
the  conditions  under  which  an  initial  crack  can  nucleate  additional  cracking  are  derived.  Only  the 
basic  equations  needed  to  calculate  the  increase  in  damage  at  the  shock  front  will  be  repeated 
here.  In  3D  the  initial  damage  is  defined  as 

Dq  =|^-(acosj)3AK,  (1) 
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where  a  is  the  radius  of  penny-shaped  cracks,  %  is  angle  describing  the  orientation  of  the  cracks 
(see  Figure.  2)  and  Ny  is  the  number  of  cracks  per  unit  volume.  In  response  to  loading,  wing 
cracks  of  length  /  grow  at  opposite  edges  of  the  initial  crack  thereby  increasing  the  damage  to 


Z)  = -j  #(/ +  acos%f  Ny. 
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(2) 


Figure  2.  Geometry  of  a  penny-shaped  crack  of  radius  a,  which  is  extended  by  wing  cracks  of 
length  /. 

The  value  of  /  is  determined  by  letting  the  wing  cracks  grow  until  the  stress  intensity  factor  at  the 
tip  decreases  to  the  fracture  toughness  of  the  medium.  The  equation  that  must  be  solved  to 
determine  the  final  amount  of  damage  is  (Ashby  and  Sammis,  1990) 
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Here  Si  and  S3  are  the  maximum  and  minimum  normalized  principal  compressive  stresses  given 
by 


Si  = 


Si  = 


^3 


(4) 
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and  Kjc  is  the  critical  stress  intensity  for  mode  I  cracks,  a  material  property.  The  constants  in  (3) 
are 


Ci 

C2 


~(;vj/2-A 

=  <cos  xf2^- - ^72 - 

(1  +  //2)  -n 


c3=  2 

C4  =  27r(cosxf  ^sTp - l-rj2 

(jV)  -M 


(5) 


In  these  last  equations,  //  is  the  coefficient  of  static  friction  and  //is  a  correction  factor  for  the 
effective  length  of  the  crack,  typically  0.45  introduced  by  Ashby  and  Sammis  (1990)  to  bring 
their  approximate  analytical  model  into  agreement  with  numerical  simulations  in  the  limit  of 
small  /.  Also,  %  is  assumed  here  to  be  45°. 

Given  the  initial  damage  D0  and  the  principal  stresses  <j\  and  03,  equation  (3)  can  be  used  to 
calculate  the  equilibrium  state  of  damage.  This  is  a  cubic  equation  in  the  damage,  which  we 
solve  numerically.  There  are  three  possible  outcomes.  At  low  stresses,  wing  cracks  do  not 
nucleate  and  we  get  no  real  solutions.  At  intermediate  stresses  we  calculate  a  value  of  D>Dn 
until  a  maximum  is  reached  above  which  we  again  get  no  real  solutions.  Ashby  and  Sammis 
(1990)  interpret  this  maximum  as  failure  since,  for  additional  loading,  damage  increases  at 
decreasing  stress,  an  unstable  condition  leading  to  shear  localization.  We  do  not  attempt  to  model 
this  post-failure  regime  beyond  identifying  it  as  the  granulated  region  described  in  the  previous 
section,  and  as  a  possible  site  for  further  weakening  by  acoustic  fluidization.  Figure  3  shows  a 
conceptual  model  of  how  damage  is  generated  in  the  stress  field  at  the  shock  front  of  the 
explosion.  In  this  case  the  principal  stresses  are  the  radial  and  hoop  stress  in  the  spherical 
geometry. 
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Figure  3.  Generation  of  fracture  damage  in  the  shock  front  of  an  underground  explosion. 


3.3  The  Damage  Mechanics  of  Frozen  Rock 

Uniaxial  compressive  and  tensile  tests  at  low  temperature  have  shown  that  freezing  can  increase 
rock  strength  by  a  factor  of  4  in  porous  rock  and  by  a  factor  of  1.8  in  crystalline  rock  (Mellor, 
1973).  We  now  interpret  these  data  using  the  micromechanical  damage  model  developed  by 
Ashby  and  Sammis  (1990)  in  which  sliding  on  preexisting  cracks  in  rock  induces  additional 
fracture  damage  and  ultimate  failure.  The  effect  of  ice  in  this  model  is  to  increase  the  effective 
coefficient  of  sliding  friction  on  the  preexisting  cracks  thus  inhibiting  the  generation  of  new 
damage  and  strengthening  the  rock.  In  addition  to  strengthening,  the  damage  model  is  also  able 
to  explain  some  of  the  subtler  phenomenology  in  the  ice  data  such  as  the  differences  between 
porous  and  crystalline  rock  and  the  progressive  strengthening  observed  to  occur  as  the 
temperature  is  lowered  from  0°C  to  -150°  C. 

The  strengthening  associated  with  low  temperatures  can  be  expected  to  significantly  reduce  the 
apparent  radius  of  an  explosion  by  suppressing  crack  growth  in  the  non-linear  damage  zone,  thus 
producing  a  higher  frequency  seismic  signature  than  that  from  an  equivalent  explosion  in  rock 
above  0°C.  Also,  the  “pulse  broadening”  phenomena  observed  for  shots  in  crystalline  rock, 
which  has  been  associated  with  dynamical  weakening  of  the  damage  zone,  should  be  suppressed 
in  frozen  crystalline  rock. 
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3.3.1  The  Uniaxial  Strength  of  Frozen  Rock 


Mellor  (1973)  measured  the  uniaxial  compressive  and  tensile  strength  of  granite,  limestone,  and 
sandstone  over  a  temperature  range  from  20°  C  to  -197°  C  for  both  water  saturated  and  air-dry 
samples.  Results  from  the  compressive  tests  are  shown  in  Fig. 4  and  those  from  the  tensile  tests  in 
Fig.  5.  We  first  make  a  few  observations  about  these  results  and  offer  a  qualitative  explanation  of 
each  in  terms  of  the  damage  mechanics  model.  A  more  quantitative  model  for  the  granite  data  is 
given  in  the  next  section. 

Observation  1:  In  compression,  the  saturated  samples  for  all  three  rock-types  show  a  strong 
increase  in  strength  beginning  at  0°C  and  continuing  to  -120°  C.  The  air-dry  limestone  and 
sandstone  samples  do  not  show  a  comparable  increase  in  strength  at  low  temperatures.  However, 
both  the  air-dry  and  saturated  granite  samples  strengthen  at  approximately  the  same  rate  with 
decreasing  temperature. 

Damage  Mechanics  Explanation:  In  the  saturated  samples  for  all  three  rock-types,  frozen  water 
inhibits  sliding  on  the  preexisting  flaws,  suppressing  damage  and  strengthening  the  samples.  The 
flow  strength  of  ice  increases  as  the  temperature  falls  below  the  freezing  point  thus  increasing 
the  apparent  coefficient  of  friction  and  strengthening  the  samples.  For  the  air-dry  limestone  and 
sandstone  specimens,  there  is  not  enough  adsorbed  water  in  the  pores  to  provide  this 
strengthening.  However,  for  granite  the  preexisting  microcracks  are  much  narrower.  The 
adsorbed  water  in  the  air-dry  samples  effectively  saturates  these  cracks,  and  hence  there  is  little 
difference  between  the  low  temperature  compressive  strength  of  saturated  and  air-dry  granite. 
Note  that  the  saturated  granite  is  significantly  weaker  at  20°  C.  This  may  be  due  to  hydrolytic 
weakening  which  is  suppressed  below  0°C. 

Observation  2:  In  tension,  the  saturated  samples  again  all  show  an  increase  in  strength 
beginning  at  0°C  and  increasing  to  lower  temperatures.  However,  in  tension,  none  of  the  air-dry 
samples  showed  a  comparable  increase  in  strength.  Even  the  air-dry  granite  samples  remain 
significantly  weaker  than  the  saturated  ones  at  low  temperature. 

Damage  Mechanics  Explanation:  Failure  in  tension  is  nucleated  at  the  largest  most  dangerously 
oriented  flaw,  which  grows  unstably  to  produce  macroscopic  failure.  Failure  in  compression  is  a 
cooperative  phenomenon  involving  structures  created  by  the  interaction  of  a  myriad  of  smaller 
flaws.  The  reason  that  the  tensile  strength  in  air-dry  granite  does  not  track  that  of  the  saturated 
granite,  is  that  the  largest  crack  (the  one  responsible  for  the  tensile  failure)  is  too  large  to  be 
saturated  by  the  films  of  adsorbed  water  that  saturate  the  myriad  of  smaller  microcracks 
responsible  for  compressive  failure. 

We  now  quantify  these  observations  by  including  the  effect  of  frozen  water  in  the  Ashby  and 
Sainmis  (1990)  damage  mechanics  model  for  granite. 
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Temperature  (C) 


Figure  4.  Strength  of  granite,  limestone,  and  sandstone  in  uniaxial  compression  at  low 
temperatures  from  Mellor  (1973).  Note  that  only  in  granite  are  the  strengths  of  saturated  and  air- 
dry  samples  nearly  the  same  at  all  temperatures. 


Temperature  (C) 


Figure  5.  Strength  of  granite,  limestone,  and  sandstone  in  uniaxial  tension  at  low  temperatures 
from  Mellor  (1973).  Note  that  the  saturated  samples  for  all  three  rock-types  show  more 
strengthening  at  low  temperatures  than  do  the  air-dry  samples. 
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3.3.2  A  Damage  Mechanics  Model  for  the  Strength  of  Frozen  Granite 

Ashby  and  Sammis  (1990)  show  that  an  initial  inclined  crack  will  nucleate  tensile  wing  cracks 
(as  illustrated  in  Fig.  2)  when 


sl  =  cls2>  +  s0 


(6) 


where  5)  and  S3  are  given  in  eqn.  (4). 

1/2 

Using  the  published  value  for  the  critical  stress  intensity  factor  Kjc=  1  MPa  m  ,  Ashby  and 
Sammis  (1990)  fit  eqn.  (6)  to  data  for  the  initiation  of  microcrack  damage  in  granite  to  detennine 
/j  in  the  range  0.55  to  0.65  and  a  crack  length  2 a  near  1  mm.  The  assumption  is  that  the  cracks 
do  not  interact  during  nucleation  so  that  eqn.  (6)  is  independent  of  the  initial  crack  density  Da . 

The  failure  curve  is  sensitive  to  the  intial  crack  density,  which  (in  3D)  is  expressed  by  the  initial 
damage  D0  defined  in  eqn.  (1).  When  eqn.  (3)  is  used  to  plot  5)  as  a  function  of  D  (for  fixed 
values  of  Da,  S3  and  the  other  parameters),  5)  has  a  maximum  value  which  Ashby  and  Sammis 
(1990)  interpret  as  the  failure  stress.  For  larger  values  of  D  beyond  this  peak,  S)  decreases.  This 
weakening  rheology  leads  to  shear  localization  and  macroscopic  failure.  Fitting  the  uniaxial 
strength  of  dry  granite  at  20  °C  using  the  values  of  a  and  ju  deduced  above  from  the  nucleation 
gives  an  initial  damage  of  Da  =  0. 1. 

We  now  ask  how  much  of  an  increase  in  //  is  required  to  produce  the  compressive  strengthening 
at  low  temperatures  observed  in  Fig.  (4).  Holding  all  the  other  parameters  constant,  we  increased 
//  to  get  the  uniaxial  strength  measured  by  Mellor  (1973)  (Fig.4).  Figure  6  shows  the  value  of  // 
as  a  function  of  temperature  required  by  this  damage  mechanics  interpretation  of  Mellor’s  (1973) 
data. 


Figure  6.  Coefficient  of  friction  on  starter  cracks  required  if  the  damage  mechanics  model  is  to 
explain  the  increase  in  compressive  strength  at  low  temperatures  in  Fig.  4. 
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3.3.3  The  Apparent  Coefficient  of  Friction  of  Frozen  Cracks 

The  coefficient  of  friction  //  is  defined  by  Amonton’s  Law  that  relates  the  shear  stress  r  required 
to  initiate  sliding  between  two  surfaces  that  are  supporting  a  normal  stress  cr 


r  =  per  (7) 

In  the  case  where  the  boundary  contains  a  fluid  with  pressure  P,  friction  is  described  by  the 
apparent  friction  law 


r  =  ju(<j-P) 


(8) 


The  coefficient  of  friction  shows  the  smallest  variation  across  different  rock  types  of  any 
mechanical  parameter.  Known  as  “Byerlee’s  Law”,  p  ~  0.6  for  virtually  all  rock  surfaces 
independent  of  hardness. 

The  asperity  model  for  friction,  initially  formulated  by  Bowden  and  Tabor  (1950,  1964)  offers  a 
simple  physical  explanation  for  the  constancy  of  p.  This  model  recognizes  that  all  surfaces  are 
rough  at  some  scale,  and  that  actual  contact  occurs  at  a  limited  number  of  high  points  called 
asperities.  Under  a  normal  load  cr  the  asperities  flow  thereby  increasing  the  true  area  of  contact, 
Ac,  until  the  stress  on  the  asperities  falls  to  the  normal  yield  stress,  crv,  of  the  rock.  If  the 

apparent  area  of  the  contacting  surfaces  is  A,  we  have 

o-A  =  ayAc  (9) 


The  shear  stress  x  required  to  initiate  sliding  is  detennined  by  the  shear  yield  stress  rv  at  which 
the  asperities  fail  in  shear,  which  gives: 


rA  =  r yAc 

Equations  (9)  and  (10)  may  be  combined  to  write 

(  \ 

Ty 

T=  — cr, 

which  is  just  Amonton’s  law  where  p  =  rv  / av . 


(10) 


(11) 


To  see  why  the  ratio  zyl<Jv  is  nearly  independent  of  rock  type,  consider  the  strength  of  a  single 
asperity  with  cross-section  Ac  that  is  loaded  to  the  compressive  yield  stress  crv .  This  asperity 
will  fail  when  the  maximum  shear  stress  reaches  ry.  As  illustrated  by  the  Mohr  Circle  in  Fig.  7 
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the  normal  yield  stress  is  related  to  the  shear  yield  stress  by  ry  / ctv  =  0.5 ,  independent  of  the 
material. 


If  the  area  between  the  asperities  is  saturated  with  water,  this  will  have  no  direct  effect  on  //,  as 
long  as  P<= 0.  However,  if  this  water  freezes,  the  shear  resistance  will  increase  to 


A 


r  =  t. 


crock 


Ar 


+  T, 


yrock  Ace  A 

Assuming  that  the  normal  load  is  supported  entirely  by  the  rock  asperities  such  that 

X 


(7  =  a 


£ rock 


yrock 

we  can  express  the  ice  strengthening  by  an  apparent  coefficient  of  friction  /J  as 


r  D-  >  A,.  ,  •  r,..  A 


_  —  yrock  crock 

r1  ~  ~  A 

<J  cr„  ,  A 


yjce  cice  =fj-\- 


yice 


X 


yrock  crock 


^y rock  A Crock 


(12) 


(13) 


(14) 


Figure  7.  Stress  state  for  a  cylindrical  asperity  loaded  to  its  compressive  yield  stress  crv.  Note 
that  Ty  =  0.5crv  where  rv  is  the  shear  strength. 

Direct  observation  of  asperities  on  surfaces  of  transparent  silicates  has  shown  that  they  deform 
by  plastic  flow  during  shear  (Dieterich  and  Kilgore,  1996).  At  the  high  stresses  and  low 
homologous  temperature  (  77 Tm  where  Tm  is  the  melting  temperature)  in  the  rock  asperities,  the 
flow  stress  is  nearly  independent  of  strain-rate  or  temperature.  However,  the  flow  stress  in  ice 
near  its  melting  temperature  is  detennined  by  thermally  activated  diffusion  limited  creep 
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mechanisms  that  are  very  sensitive  to  both  temperature  and  strain-rate.  This  is  why  //  rises 
rapidly  at  temperatures  just  below  0°C  in  Fig.  6.  In  equation  (14),  this  behavior  is  caused  by  the 
rapid  increase  in  ?yice  with  falling  temperature  below  0°C.  At  very  low  temperatures,  flow  in 

the  ice  will  also  be  controlled  by  glide  and  /i  will  be  almost  independent  of  temperature  and 
strain-rate. 

Flow  in  ice  near  its  melting  point  is  more  complicated  than  high  temperature  flow  in  rock.  There 
is  a  poorly  understood  change  in  activation  energy  at  about  -10°C  thought  to  be  associated  with 
the  suppression  of  pressure  melting  at  grain  boundaries  below  this  temperature  (Goodman  et  ah, 
1981;  Duval  et  ah,  1983).  We  will  take  the  analysis  in  this  paper  one  step  further  by  assuming  a 
generic  form  for  this  high  temperature  power-law  flow  in  ice,  the  strain  rate  b  is  given  by 

N=  Bzn  Qxp(-Q/RT)  (15) 

Solving  for  the  flow  stress  gives 


^ yice 


r  ^\l/n 

By 


exP  (Q/nRT) 


Using  this  in  eqn.  (14)  and  combining  constants  gives 

/J=  ju+  Dcxp(Q/ riRTj 

or,  equivalently, 


In (ju'-ju)  =  In  D  +  Q/nRT 


(16) 


(17) 

(18) 


T  (C) 


227  -23  -106  -148  -173  -190  -202 


1/T(K) 


Figure  8.  Plot  of  equation  (17)  in  the  text.  The  slope  in  the  thermally  activate  region  gives 
Q/nR  =  714 .  At  temperatures  below  about  1 10°C  the  mechanisms  changes  to  glide  controlled 
plasticity. 
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In  Fig. 8,  In (//’-//)  is  plotted  as  a  function  of  1  IT.  The  linear  portion  has  a  slope  of  Q/nR  =  714. 
Below  8.3°C,  the  activation  energy  for  power  law  creep  is  Q  ~  80  KJ  mole' ^Goodman  et  ah, 
1981).  This  implies  a  creep  power  of 


8xl04 

"'(714X8.3) 


13.5 


(19) 


This  is  significantly  higher  than  n  ~  3  observed  at  stresses  below  IMPa.  However,  n  for  ice  is 
observed  to  rise  rapidly  at  stresses  above  1  MPa  (Goodman  et  ah,  1981),  so  this  may  simply 
reflect  a  higher  flow  stress  on  the  asperities. 

3.3.4  The  Effect  of  Ice  on  the  Elastic  Properties  of  Frozen  Rock 

The  effect  of  ice  in  the  inclined  crack  in  Fig.  2  is  to  either  totally  immobilize  it,  thus  reducing 
Ny  in  eqn.  (1),  or,  if  saturation  is  not  total,  to  form  ice  bridges  thus  reducing  a  in  eqn.  (1)  while 
increasing  Ny .  The  net  effect  of  both  is  to  reduce  the  initial  damage  D0. 


The  elastic  moduli  of  rock  are  extremely  sensitive  to  D0.  Figure  9,  from  O’Connell  and 
Budiansky  (1974)  shows  the  effect  of  changing  D0  on  the  P  and  S  wave  velocities  in  dry  and 

water  saturated  rock.  Note  that  the  x-axis  in  Fig. 9  is  £  =  n[c£>\  where  N  is  the  number  of  cracks 

per  unit  volume  and  a  is  the  half-length  of  the  crack.  Comparison  with  eqn.  (1)  shows  that  s  can 
be  written  in  terms  of  the  initial  damage  as 


£  = 


4  7t 


1  33 


vcos  XJ 


Da  «  0.68Do 


(20) 


The  effect  of  freezing  water  in  the  cracks  is  to  move  to  the  left  (toward  lower  damage)  on  the 
curves  in  Fig. 9.  This  will  produce  an  increase  in  elastic  wave  velocity.  The  effect  is  larger  for  S 
waves  than  for  P  waves  in  saturated  rock. 


Figure  9.  The  effect  of  fractures  on  the 
elastic  wave  velocities  in  wet 
and  dry  rock  (from  O’Connell 
and  Budiansky,  1974).  The 
crack  density  parameter  £  is 
closely  related  to  the  damage 
parameter  in  the  Ashby  and 
Sammis  (1990)  damage 
mechanics,  differing  only  by  a 
constant  (see  eqn.20). 
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3.3.5  The  Effect  of  Strain-rate  on  the  Strength  of  Frozen  Rock 

The  effect  of  strain-rate  on  the  strength  of  frozen  rock  is  illustrated  by  the  deformation- 
mechanism  map  for  ice  in  Fig.  10.  Temperatures  near  melting  are  at  the  right  edge  of  the  map 
where  it  can  be  seen  that  increasing  the  strain  rate  produces  a  significant  increase  in  the  flow 
stress.  At  an  effective  shear  stress  of  about  10  MPa  this  trend  is  reversed  by  the  onset  of  pressure 
melting  which  weakens  the  ice. 


temperature  °c 
-200  -100 


-  Onset  of  pressure  melting 

Increase  in  shear  strength  with 
Increasing  strain-rate. 


Figure  10.  Deformation-mechanism  map  for  ice.  Note  that  increasing  the  strain-rate  by  5  orders 
of  magnitude  near  OC  increases  the  shear  strength  (effective  shear  stress)  by  a  factor  of  10.  This 
increase  in  strength  is  limited  by  pressure  melting  that  occurs  at  an  effective  shear  stress  of  about 
10  MPa.  (From  Frost  and  Ashby,  1982). 

3.4  The  Effect  of  Ice  in  Rock  on  Seismic  Radiation  from  an  Explosion 

In  order  to  model  the  effect  of  ice  on  the  damage  generated  by  an  explosion  and  on  the  seismic 
radiation,  we  require  the  principal  stresses  generated  by  the  explosion  as  a  function  of  distance 
and  time.  This  is  made  difficult  by  the  existence  of  the  nonlinear  processes  between  the  cavity 
radius  and  the  effective  elastic  radius,  beyond  which  the  assumptions  of  ordinary  linear  elasticity 
are  valid.  Sophisticated  computer  codes  have  been  developed  which  include  hydrodynamic 
effects,  shock  waves,  and  nonlinear  equations  of  state  (see,  for  example,  Rodean,  1971;  King  et 
ah,  1989;  Glenn,  1993;  Glenn  and  Goldstein,  1994,  for  discussion  and  further  references).  We 
use  here  an  approximate  method  to  calculate  the  stresses  surrounding  an  explosion  that  is  based 
on  the  equivalent  elastic  method  developed  by  earthquake  engineers  to  model  the  nonlinear 
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behavior  of  soils  that  that  occurs  during  strong  ground  motion.  The  central  idea  is  to  make  the 
material  properties  a  function  of  the  stress  in  the  outward  propagating  pressure  pulse  and  then  to 
adjust  these  material  properties  in  an  iterative  process  until  the  appropriate  values  are  present  at 
all  distances  from  the  source.  In  effect,  the  nonlinear  stress-strain  behavior  is  approximated  by  a 
series  of  linear  relationships  that  change  with  the  level  of  stress.  The  present  formulation, 
described  by  Johnson  (1993),  relates  density  and  bulk  elastic  properties  to  the  peak  pressure  and 
shear  and  anelastic  properties  to  the  maximum  shear  strain. 

The  details  of  this  model  are  published  in  Johnson  and  Sammis  (2001)  and  will  not  be  repeated 
here. 

In  Johnson  and  Sammis  (2001)  we  modeled  the  1  kt  chemical  explosion  detonated  in  September 
1993  as  part  of  the  Non-Proliferation  Experiment  (NPE)  (see  Denny,  1994)  and  found  that  the 
equivalent  elastic  model  gave  a  good  representation  of  the  near  field  amplitudes.  We  now  ask 
how  these  results  would  change  if  we  increased  the  P  and  S  wave  velocities,  the  coefficient  of 
friction,  and  the  strength  to  simulate  the  mechanical  properties  if  the  cracks  were  filled  with  ice. 

Figure  1 1  shows  the  effect  on  the  radial  distribution  of  fracture  damage  in  the  non-linear  source 
region  when  the  coefficient  of  friction  is  increased  from  //=  0.6  to  ju=  1 .0.  For  the  higher 
coefficient  of  friction,  damage  is  suppressed  close  in  and  enhanced  further  out.  This  is  because 
the  stresses  do  not  fall  off  as  fast  with  distance  due  to  the  smaller  damage  close  in.  Radial  cracks 
are  thus  expected  to  extend  out  an  additional  10  meters  in  the  frozen  rock  (about  an  additional 
20%). 
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Figure  11.  The  effect  on  the  radial  distribution  of  damage  when  the  coefficient  of  static  friction 
is  increased  from  /u  =  0.6  to  /u  =  1 .0  in  the  NPE  simulation  by  Johnson  and  Sammis 
(2001). 
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Figure  12  shows  the  effect  of  increasing  the  compressive  strength  on  the  radial  distribution  of 
damage.  In  the  equivalent  elastic  source  model  the  compressive  strength  is  expressed  as  a 
reference  strain  defined  as 


£y 


^max 

Mo 


(21) 


In  Fig.  12,  the  reference  strain  is  doubled  (corresponding  to  ice  in  granite)  and  multiplied  by  4 
(corresponding  to  ice  in  limestone).  The  effect  is  very  similar  to  that  of  increasing  the  coefficient 
of  static  friction  -  there  is  more  damage  and  it  extends  an  additional  10  meters. 
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Figure  12.  The  effect  on  the  radial  distribution  of  damage  when  the  non-linear  reference  strain 
is  increased  by  factors  of  2  and  4  in  the  NPE  simulation  by  Johnson  and  Sammis 
(2001).  This  is  equivalent  to  increasing  the  compressive  strength  by  a  factor  of  2 
(simulating  ice  in  granite)  and  4  (simulating  ice  in  limestone). 


Figure  13  summarizes  the  effect  of  ice  in  the  cracks  on  the  amplitude  of  seismic  radiation  in  the 
far-held.  As  argued  above,  ice  is  expected  to  increase  the  initial  wave  velocities  and  the  strength 
of  the  source  rock.  These  ice-induced  changes  all  decrease  both  the  amplitude  and  reduced 
velocity  potential  in  the  linear  elastic  regime. 
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Figure  13  The  effect  of  ice  in  the  cracks  on  the  reduced  velocity  potential  and  far- Held 

displacement.  The  solid  circle  is  the  velocity  potential  and  displacement  at  a  radial 
distance  of  650  m  from  the  NPE  explosion  calculated  by  Johnson  and  Sammis  (2001). 
The  other  symbols  show  the  effect  of  changing  either  the  initial  velocity  or  strength  of 
the  rock.  Note  that  effect  of  increasing  either  the  initial  seismic  velocity  or  the  strength 
is  to  reduce  the  amplitude  of  the  seismic  radiation  in  the  far  field  making  the  explosion 
appear  smaller  than  an  equivalent  explosion  in  rock  above  the  freezing  temperature  of 
water. 

3.5  Field  Measurements  of  the  Effect  of  Ice  on  Seismic  Coupling:  The  Alaska 
Experiment 

In  collaboration  with  Weston  Geophysical  and  New  England  Research  a  series  of  identical 
chemical  explosions  were  detonated  above  and  below  the  permafrost  layer  in  crystalline  bedrock 
at  a  site  near  Fairbanks  Alaska  (Fig.  14).  The  details  of  this  experiment  are  given  in  Bonner  et  al. 
[2007].  Results  summarized  here  include  source  spectra  for  several  explosions  at  both  sites  and 
the  seismic  P  and  S  velocities  at  each.  Ultrasonic  P  and  S  wave  velocities  in  laboratory  samples 
from  each  test  site  are  presented  here  to  supplement  these  field  data.  The  primary  result,  shown 
in  Figure  15,  was  that  the  spectral  levels  of  radiation  from  both  sites  were  comparable,  but  the 
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comer  frequency  of  radiation  from  the  frozen  rock  site  was  higher.  This  higher  corner  frequency 
would  lead  to  the  erroneous  conclusion  that  the  explosion  in  frozen-rock  was  smaller. 

The  two  test  sites  are  shown  on  the  map  (Fig.  14a)  and  cross-section  (Fig.  14b).  The  first, 
hereafter  referred  to  as  the  “dry  site”,  was  located  above  the  water  table  in  a  mica  schist.  The 
second,  hereafter  referred  to  as  the  “frozen  site”,  was  located  below  the  permafrost  in  a  gneiss. 
Details  of  the  seismic  instrumentation  and  deployment  are  given  in  Bonner  et  al.  [2007]. 
Measurement  of  the  seismic  velocities  at  the  dry  site  gave  vP=  2,600  m/s  and  v$  =  1,600  m/s, 
while  measurements  at  the  frozen  site  gave  vp  =  3,600  m/s  and  =  2,200  m/s.  The  yield  and 
depth  of  burial  of  each  shot  are  summarized  in  Table  1. 

Table  1.  Location,  yield,  and  depth  of  burial  of  explosions  near  Fairbanks  Alaska 
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Figure  14.  Upper  panel 
shows  the  location  of  frozen 
and  unfrozen  test  sites. 
Lower  panel  shows  cross- 
section  A- A’  as  indicated  on 
the  upper  map. 
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Near-Source  Spectra 


Figure  15.  Displacement  spectra  for  explosions  in  unfrozen  dry  schist  (UN)  and  frozen  gneiss 
(FR).  Note  the  spectral  peak  and  higher  corner  frequency  for  shots  at  the  frozen  site.  (Figure 
courtesy  of  Jessie  Boner  at  Weston  Geophysical). 


In  order  to  make  a  quantitative  interpretation  of  this  result,  we  wish  to  model  the  Alaska 
explosions  using  the  Johnson  and  Sammis  (2001)  source  model.  This  model  incorporates  the 
Ashby  and  Sammis  (1990)  micromechanical  damage  mechanics,  which  as  shown  above,  gives  a 
good  description  experimental  data  for  the  temperature  dependence  of  strength  in  frozen 
crystalline  rock  (Sammis  and  Biegel,  2004,  2005).  The  basic  idea  is  that  ice  in  the  cracks 
strengthens  the  rock  and  increases  its  elastic  stiffness  by  reducing  the  initial  damage  Da.  In  our 
source  model,  this  results  in  less  explosion-induced  damage  in  the  frozen  rock,  and  a  higher 
comer  frequency  -  as  observed. 

We  use  the  theory  developed  by  O’Connell  and  Budiansky  (1974)  to  estimate  D0  from  the 
seismic  velocities  measured  by  Weston  Geophysical  at  each  field  site.  Figure  9  above  shows  the 
ratios  vp  /vp  and  /v$  as  functions  of  a  crack  density  parameter,  which  they  define  as 

s  =  Ny(a)  and  which  is  simply  related  to  D0  in  equation  (20)  .The  quantities  vP  and  vs  are 
velocities  in  the  fractured  rock,  while  vP  and  vs  are  velocities  in  the  unfractured  rock. 

In  order  to  use  Figure  9  to  find  the  crack  density  parameter  £’(and  hence  the  initial  damage  D0) 
from  the  seismic  velocities  ( vp  and  v$)  at  each  test  site,  we  need  to  know  the  un- fractured 
velocities  ( vp  and  v^)  for  the  source  rock  at  each  site  (the  schist  and  the  gneiss).  These  we 
estimate  using  the  ultrasonic  velocities  measured  in  laboratory  samples  of  each  rock  type  by  New 
England  Research  (NER). 
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For  the  schist  source  rock  from  the  dry  (un-frozen)  site,  NER  measured  P  and  S  wave  velocities 
under  dry,  saturated,  and  frozen  conditions,  at  confining  pressures  ranging  from  5  to  70  MPa. 
Table  2  summarizes  the  values  at  5  and  at  70  MPa  -  values  at  intermediate  pressures  vary 
monotonically  between  these  limits. 


Table  2.  Ultrasonic  data  for  schist  source  rock  from  the  dry  unfrozen  site. 


P  (MPa) 

Vp  (m/s) 

V5  (m/s) 

V 

Vp  (m/s) 

Vs*  (m/s) 

Dry 

5 

4150 

2800 

0.082 

5253 

3294 

70 

4750 

3200 

0.085 

5723 

3556 

Saturated 

5 

5200 

3000 

0.25 

5306 

3261 

70 

5450 

3400 

0.18 

5561 

3579 

Frozen 

(-8C) 

5 

5800 

3700 

0.16 

70 

5700 

3700 

0.14 

*Unfractured  reference  velocities  vP  and  vs  found  using  <5pand  s$  from  Table  2. 


The 

that 


last  two  columns  of  Table  2  contain  estimates  of  the  un-fractured  P  and  S  wave  velocities 

were  made  using  the  calculations  summarized  in  Table  2.  In  Table  2,  values  of  ^(satd  anci 

vp(clry) 


ES'CraT)  arc  calculated  from  the  data  in  Table  1 .  Note  that  A  p  an^  that 

vs(dry)  vp(dry)  vP(dry)/vP 

_  Es'CraT)  A  $  anc^  ratios  on  the  right  hand  side  of  these  equations  can  be  found 

vs (dry)  vs(dry)/vs 

as  a  function  of  s'm  Figure  9  and  are  plotted  explicitly  in  Figure  16.  The  values  of  s'm  Table  2 
were  found  using  Figure  16  as  indicated  on  the  figure.  Once  s  was  determined,  the  ratios  v  /v 
were  found  from  Figure  9  and  appear  in  the  last  four  columns  of  Table  2.  These  ratios  were  used 
together  with  the  measured  seismic  velocities  to  calculate  the  un-fractured  velocities  in  the 
gneiss  samples,  which  are  give  in  the  last  two  columns  of  Table  2.  Note  that  the  un-fractured 
velocities  found  by  this  method  are  very  similar  to  the  velocities  measure  in  the  frozen  samples. 
This  observation  supports  our  assumption  that  ice  immobilizes  the  fractures. 


Table  3.  Ultrasonic  velocity  ratios  used  to  find  fracture  parameter  s  in  gneiss  lab  samples. 


P  (MPa) 

vP(sat.) 

vp(dry) 

vs(sat.) 
vs  (dry) 

£P 

% 

Vp 

Vp 

(dry) 

Vp 

Vp 

(sat.) 

US 

VS 

(dry) 

US 

(sat.) 

5 

1.25 

1.07 

0.18 

0.16 

0.79 

0.98 

0.85 

0.92 

70 

1.15 

1.06 

0.13 

0.15 

0.83 

0.98 

0.90 

0.95 
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Figure  16.  The  ratio  of  saturated  to  dry  velocities  as  a  function  of  the  crack  parameter  a.  The 
curves  were  calculated  from  Figure  2.  The  points  are  ratios  of  the  ultrasonic  data  for  schist  from 
Table  2.  Note  the  decrease  in  crack  parameter  as  confining  pressure  closes  cracks. 


For  the  gneiss  source  rock  from  the  frozen  site,  NER  measured  P  and  S  wave  velocities  under 
saturated  and  frozen  conditions,  also  at  confining  pressures  ranging  from  5  to  70  MPa.  Table  4 
summarizes  the  values  at  5  and  at  70  MPa  -  values  at  intennediate  pressures  vary  monotonically 
between  these  limits. 


Table  4.  Ultrasonic  velocities  for  gneiss  source  rock  from  the  frozen  site. 


P  (MPa) 

Vp  (m/s) 

V5  (m/s) 

V 

Saturated 

5 

5450 

3050 

0.27 

70 

5550 

3200 

0.25 

Frozen 

(-8C) 

5 

5800 

3400 

0.24 

70 

5800 

3450 

0.23 
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Since  no  measurements  were  made  under  dry  conditions,  we  can  not  use  the  same  method  we 
used  to  estimate  the  un-fractured  velocities  in  the  schist.  However,  based  on  the  results  for  the 
schist,  we  can  assume  that  the  frozen  velocities  are  a  reasonable  approximation  to  the  un¬ 
fractured  velocities. 


Having  estimated  the  un-fractured  velocities  vp  and  v$  for  both  source  rocks,  it  is  now  possible 
to  estimated  the  fracture  parameter  at  each  test  site.  The  seismic  field  velocities  and  implied 
fracture  parameters  at  each  site  are  summarized  in  Table  5. 


Table  5.  Ultrasonic  velocities  for  schist  source  rock  from  the  dry  unfrozen  site. 


Vp  (m/s) 

(m/s) 

V 

Vp_ 

Vp 

Us 

£P 

% 

(e) 

A) 

Unfrozed  Site 

2700 

1700 

0.17 

0.47 

0.46 

0.41 

0.45 

0.43 

1.80 

Frozen  Site 

3800 

2400 

0.17 

0.66 

0.65 

0.28 

0.35 

0.31 

1.30 

Note  that  for  both  P  and  S  waves,  the  fracture  parameter  s,  and  the  equivalent  initial  damage  Da , 
is  larger  at  the  unfrozen  site,  consistent  with  the  hypothesis  that  ice  reduces  the  fracture  density 
in  the  source  rock.  When  these  values  of  Da  are  used  in  the  Johnson  and  Sammis  (2001)  source 
model,  the  unfrozen  site  has  a  larger  elastic  radius  and  lower  corner  frequency  than  the  frozen 
site,  as  observed  in  the  Alaska  frozen  rock  experiment. 

4.  RESULTS  AND  DISCUSSION 

Laboratory  studies  of  the  effects  of  ice  on  the  material  properties  of  rock  find  that  it  increases 
the  elastic  stiffness  and  strength.  These  effects  are  observed  to  be  both  temperature  and  strain- 
rate  dependent,  with  the  stiffness  and  strength  increasing  at  lower  temperatures  and  higher  strain- 
rates.  These  observations  can  be  explained  using  a  micromechanical  damage  mechanics  where 
the  primary  effect  of  ice  in  the  cracks  is  to  decrease  the  crack  density  and  increase  the  coefficient 
of  sliding  friction.  When  applied  to  a  field  experiment  in  which  identical  explosions  were 
detonated  in  frozen  and  unfrozen  rock,  the  suppression  of  damage  by  ice  in  the  cracks  explains 
the  observation  of  a  smaller  elastic  radius  at  the  frozen  site. 

5.  CONCLUSIONS 

Frozen  rock  can  have  a  significant  effect  on  seismic  coupling  of  underground  nuclear 
explosions  detonated  in  the  permafrost  layer.  By  suppressing  fracturing  in  the  non-linear  source 
region,  such  explosions  have  a  smaller  elastic  radius  than  those  detonated  in  unfrozen  rock.  They 
radiate  more  high  frequency  energy  and  appear  to  have  a  lower  yield. 
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